Towards a Stable Numerical Evolution of Strongly Gravitating Systems in General 

Relativity: The Conformal Treatments 
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We study the stability of three-dimensional numerical evolutions of the Einstein equations, com- 
paring the standard ADM formulation to variations on a family of formulations that separate out 
the conformal and traceless parts of the system. We develop an implementation of the conformal- 
traceless (CT) approach that has improved stability properties in evolving weak and strong gravi- 
tational fields, and for both vacuum and spacetimes with active coupling to matter sources. Cases 
studied include weak and strong gravitational wave packets, black holes, boson stars and neutron 
stars. We show under what conditions the CT approach gives better results in 3D numerical evo- 
lutions compared to the ADM formulation. In particular, we show that our implementation of the 
CT approach gives more long term stable evolutions than ADM in all the cases studied, but is less 
accurate in the short term for the range of resolutions used in our 3D simulations. 
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I. INTRODUCTION 

Three dimensional (3D) numerical relativity is an im- 
portant technique for exploring the strong field dynam- 
ics in realistic astrophysical phenomena involving black 
holes and neutron stars. It is expected to play a role in 
analyzing gravitational waveforms to be observed soon, 
one expects, with the new generation of gravitational 
wave detectors going online worldwide in the next few 
years 0,^]. However, progress in 3D numerical relativ- 
ity, which has traditionally been based on the Arnowitt- 
Deser-Misner (ADM) || system of evolution equations, 
has been slow. This is not only because of the immense 
computational difficulties that 3D simulations represent, 
but to a large extent it is due to severe instabilities often 
encountered during such simulations. Presently there is 
no complete understanding of the causes of these insta- 
bilities in numerical evolutions of the ADM equations. 
This has prompted much recent effort in developing al- 
ternative formulations of the 3+1 Einstein equations. 

In this and a companion paper 0] we focus on an alter- 
native approach based on a conformal decomposition of 
the metric and the trace-free components of the extrin- 
sic curvature. The conformal-tracefree (CT) approach 
was first devised by Nakamura in the 1980's in 3D cal- 
culations [|5|,[|, and then modified and applied to work 
on gravitational waves 0, and on neutron stars 
This approach was not taken up by others in the commu- 
nity until a recent paper by Baumgarte and Shapiro [fiof , 



where a similar formulation was compared with the stan- 
dard ADM approach and shown to be superior, in terms 
of both accuracy and stability, on tests involving weak 
gravitational waves, with geodesic and harmonic slic- 
ing. In a followup paper, Baumgarte, Hughes, and 
Shapiro fi"lf| applied the same formulation to systems 
with given (analytically prescribed) matter sources, and 
found similar stability properties. More recently fully 
hydrodynamical simulations employing the CT approach 
have been reported in ]T^-|l^| in the context of collapse of 
rapidly-rotating (isolated) neutron stars and coalescence 
and merger of binary neutron stars. As we were prepar- 
ing this manuscript we have also become aware of work 
by Lehner, Huq and Garrison Jl5[ where a comparison 
between the ADM and CT formulation in spherical sym- 
metry has been carried out in the context of black hole 
excision. 

In the companion paper [Q we perform an analytic in- 
vestigation of the stability properties of the ADM and the 
CT evolution equations. Using a linearized plane wave 
analysis, we identify features of the equations that we 
believe are responsible for the difference in their stability 
properties. 

In this paper we report the results of simulations of 
weak and strong gravitational wave packets, black holes, 
boson stars and neutron stars in various slicing condi- 
tions, including maximal slicing and a family of algebraic 
slicings, and compare the results obtained by the ADM 
and CT equations in different implementations. We be- 
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gin with a brief presentation of the relevant equations in 
Sec. ||. We then discuss the results of our numerical sim- 
ulations in_scction III . We consider vacuum spacetimes in 
section III A , and matter spacetimes in section [II B . In 
section III A 1, we describe the various implementations 



of the CT equations using gravitational wave spacetimes 
as an example. We identify two particular implementa- 
tions, which we call AFA and AF2, that give the best per- 
formance in long term evolutions. The essence of these 
implementations, is to "actively force" (AF, see below) 
the trace of the conformally rescaled extrinsic curvature 
(AFA), and for maximal slicing also the trace of the ex- 
trinsic curvature (AF2), to zero in each step of the nu- 
merical evolution. In the sections that follow, we focus 
on comparing the AFA and AF2 implementations to the 
results of the ADM equations for evolutions of strong 
field systems including black holes, boson stars and neu- 
tron stars. We demonstrate that for this wide range of 
systems, these two implementations of the CT equations 
always lead to more stable long term evolutions. How- 
ever, we also find that for a given resolution, the ADM 
results are often more accurate than the CT results at 
early times, before the instabilities become apparent. We 
conclude with section IV . A study of the stability proper- 
ties of the iterative Crank-Nicholson (ICN) scheme, used 
for the spacetime evolution of the simulations presented 
in this paper, can be found in the Appendix. 



II. FORMULATION 

We start reviewing briefly the formulations used for 
the comparisons. 

The standard ADM equations are p6| : 



d 

— 7^ = -2aKy, 
d 

dt J 



(1) 



DiDja + a \ Rij + KKij 
2K ik K k j - ^Rij) , (2) 



The conformal, trace-free reformulations of these equa- 
tions make use of a conformal decomposition of the three- 
metric, and the trace-free part of the extrinsic curvature. 
Here we follow closely the presentation of Ref . jnj . The 
conformal three-metric %j is written as 



_ -40 



with the conformal factor chosen to be 



S ^= 7 1 / 3 ^det( 7y -) 1 



/3 



(5) 



(6) 



In this way the determinant of 7y is unity. The trace-free 
part of the extrinsic curvature K%j , defined by 



A%j — ^ij 



(7) 



where K = is the trace of the extrinsic curvature, 

is also conformally decomposed: 



A* 



(8) 



So far, these are just definitions of new variables, with 
no clear motivation for their introduction. Evolution 
equations for these new quantities are easy to find, and 
we summarize here the Baumgarte-Shapiro [ [To[ discus- 
sion on these equations, but with an emphasis on the 
possible numerical implications of various choices one can 
make. 

The evolution equations for the conformal three- 
metric 7y , and its related conformal factor <fi are trivially 
written as 
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The evolution equation for the trace of the extrinsic 
curvature K can easily be found to be 



dt ' 3 



AnA ij + -K 



\{P + S) 



with 



(11) 



d_ 

dt 



dt-C 



(3) 



and where Cp is the Lie derivative with respect to the 
shift vector Here Rij is the Ricci tensor and Di the co- 
variant derivative associated with the three-dimensional 
metric jij. The 4-dimensional Ricci tensor ^'Rij is usu- 
ally written in terms of the energy density p and stress 
tensor SV,- of the matter as seen by the normal (Eulerian) 
observers: 



where the Hamiltonian constraint was used to eliminate 
the Ricci scalar. 

For the evolution equation of the trace-free extrinsic 
curvature A^ there are many possibilities. A trivial ma- 
nipulation of Eq. @ yields: 



dt ■ 



-4<j> 



-DiDjct 



+a(KA ii -2A u A l 



(12) 
(13) 



(4)d.. _ 



8tt 



(4) 



but as shown previously [j7 10 there are many ways to 
write several of the terms, especially those involving the 
Ricci tensor. For example, one could eliminate the Ricci 
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scalar R again through the use of the Hamiltonian con- 
straint. 

With the conformal decomposition of the three-metric, 
the Ricci tensor now has two pieces, which we write as 

Rn =Ri3+R%- (14) 

The "conformal-factor" part Rfj is given directly by 
straightforward computation of derivatives of 4>: 

Rf 3 = -2DiDj(l) ~ 2j l] D l D l cj) (15) 
+±b i <t>D j <t>-4r l i j D l <t> Di4>, (16) 

while the "conformal" part Rij can be computed in the 
standard way from the conformal three-metric jij. To 
simplify notation, it is convenient to define what Ref. |l(]] 
calls the "conformal connection functions" : 

r-.= ^ k r jk = -f j d , (i7) 

where the last equality holds if the determinant of the 
conformal three-metric 7 is actually unity (notice that 
this should be true analytically, but may not be numeri- 
cally). 

Using the conformal connection function, the Ricci ten- 
sor can be written: 

Rij = --jl lm iij,lra + lk(idj)t k + f k T (ifi k 

+l lm ( 2r f(*f j)km + f imf feij) ■ ( 18 ) 

Here again, one has choices in how the terms involv- 
ing the conformal connection functions T l are computed. 
A straightforward computation based on the Christof- 
fel symbols could be used (and usually is in standard 
ADM formulations), but this approach leads to deriva- 
tives of the three-metric in no particular elliptic form. 
One would like to see an elliptic form as the principal 
part of this expression, as it brings the 7^ — Ay sys- 
tem a step closer to being hyperbolic. Thanks to the 
definition of the r l 's, an explicitly elliptic operator is 
singled out. However, if the terms involving the T z are 
evaluated directly in terms of derivatives of the three- 
metric, this elliptic operator serves no special purpose, 
as other second derivatives appear through derivatives of 
the T l which spoils the elliptic nature of the operator as 
a whole. If, on the other hand, the T l are promoted to 
independent variables, for which evolution equations can 
be derived, then the expression for the Ricci tensor re- 
tains its elliptic character. The price to pay is that one 
must now evolve three additional quantities in the evolu- 
tion system. Whether this has any numerical advantage 
will depend on details of the implementation, and will be 
discussed below. 

Following this argument of promoting the T z to inde- 
pendent variables, it is straightforward to derive their 
evolution equation: 



dt ~ dxA ZaA Zl p > m 

+ r s f 3 (3\ l +P l l J A ). (19) 

However, again there is a choice one can make in writ- 
ing this evolution equation; as pointed out in Ref. ll(J 
it turns out that the above choice leads to an unstable 
system. A choice which will be shown to be better can 
be obtained by eliminating the divergence of A 13 with the 
help of the momentum constraint: 

J^f 1 = -2A ij a d + 2a{f) k A k 3 

-\rf j K tj - f j Sj + (i.r-'o ; ) 

-^j - 2f nU P% + \i l3 P\) ■ (20) 

With this reformulation, in addition to the evolu- 
tion equations for the conformal three-metric 7y (|^) 
and the conformal-traceless extrinsic curvature variables 
Aij (|l^) , there are evolution equations for the conformal 
factor (j) (|To|), and the trace of the extrinsic curvature 
K (pi]). If the f l are promoted to the status of fun- 
damental variables, as in Ref. @, they can be evolved 
with (p0|). (Note that the mixed first and second order 
evolution system for {</>, K, 7^, A^, T 1 } is not in any im- 
mediate sense hyperbolic ]l7j]). In the original formula- 
tion of Shibata and Nakamura H , the auxiliary variables 
Fi = — Y^j are used instead of the T l , and the final 
system of equations is somewhat more complicated. 

Ref. |l8| shows that the CT system can also be in- 
terpreted as a "conformal second-order" version of the 
Bona-Masso system with 2Vj = — (Tj + 88^). 

A. Gauge 

Systems of the CT type have been investigated with 
various slicing conditions in the past. The paper of 
Baumgarte and Shapiro considered geodesic and har- 
monic slicing, while earlier work by Shibata and Naka- 
mura, and the more recent paper by Baumgarte, Hughes, 
and Shapiro |llj have also considered maximal slicing. 
Here we have studied maximal slicing and a number of 
algebraic slicings, and used them with different imple- 
mentations of the CT equations, on numerical evolutions 
of many different spacetimcs. 

Maximal slicing has the property that K = 0, leading 
to an elliptic equation for the lapse 

V 2 a = a [KijK ij + 4tt (p + S)] . (21) 

Notice that in maximal slicing the evolution equations 
for <fi and K become simply 

d(j)/dt = 0, dK/dt = 0. (22) 
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The algebraic slicings that we will consider here cor- 
respond to the family originally introduced by Bona and 
Masso , building on earlier work of Bernstein PQ| 

Ct a = -f(a)a 2 K, (23) 

with f(a) > but otherwise arbitrary. This family con- 
tains many well known slicing conditions. For example, 
taking / = 1 we recover the "harmonic" slicing condition, 
which after a trivial integration becomes 

a = F{x l )+ 1 1 ' 2 , (24) 

with F an arbitrary function of space. The name "har- 
monic" slicing comes from the fact that it corresponds to 
the choice of a harmonic time coordinate 

at = 0. (25) 

Another useful slicing condition is obtained by taking 
/ = N/a. This corresponds to the generalized "1+log" 
slicing condition Jig ] which after integration becomes 

a = F{x l )+logj N/2 . (26) 

(There is in fact some inconsistency in terminology as to 
whether the N = 1 or the N = 2 case corresponds to the 
standard "1+log" slicing; different choices being made by 
different authors.) 

These type of algebraic slicings have an advantage over 
maximal slicing in terms of computational efficiency: ft 
is much faster to integrate an evolution equation for the 
lapse than to solve an elliptic equation. On the other 
hand, such algebraic slicings are prone to the develop- 
ment of gauge pathologies (|^J2^] . The possibility of the 
appearance of such pathologies when using algebraic slic- 
ings should always be kept in mind, as a gauge pathology 
can easily be confused with a numerical instability: one 
can lose a lot of sleep trying to cure an "instability" that 
is in fact a true solution of our system of differential equa- 
tions. 

To finish discussing our choice of gauge, we need to 
mention the fact that all the simulations described here 
have been carried out with the shift vector set to zero. 



allows waves coming from any arbitrary direction to leave 
the grid with no reflections, and second, there does not 
even exist a clear way to define what a wave is in gen- 
eral relativity except at asymptotic infinity. In practice, 
what one looks for is a condition that remains stable and 
allows some "wave-like" solutions to leave the grid with- 
out introducing large reflections at the boundaries. The 
amount of artificial reflection that results typically de- 
pends on the specific form of the boundary condition, 
and on the direction of motion of the wave fronts as they 
hit the boundary |23| ]. 

Since in this paper we are interested in the question of 
the stability of the interior evolution, we will not worry 
too much about the boundary conditions, and we will 
limit ourselves to describing a few conditions that we 
have found to work well in practice. The conditions we 
have used are the following: 

• Static boundary condition: The evolved variables 
are simply not updated at the boundary, and re- 
main with their initial values there. This condition 
is very bad at handling waves since it reflects every- 
thing back in, but it can be very useful when study- 
ing situations that are supposed to remain static (as 
are some of the systems studied below) , and where 
all the dynamics comes from numerical truncation 
errors. 

• Zero-order extrapolation or "flat" boundary con- 
dition: After evolving the interior, the value of a 
given variable at the boundary is simply copied 
from its value one grid point in (along the nor- 
mal direction to the boundary). This condition 
allows for some dynamics at the boundaries, and 
is somewhat better at absorbing waves than the 
static boundaries, but it still introduces a consid- 
erable amount of reflections. 

• Sommerfeld or "radiative" boundary condition: In 
this case we assume that the dynamical variables 
behave like a constant plus an outgoing radial wave 
at the boundaries, that is: 

f(x\t) =f + u (r-t)/r, (27) 



B. Boundary conditions 

In standard 3+1 numerical simulations, the computa- 
tional domain covers only a finite region of space. One 
must therefore apply some sort of artificial boundary con- 
dition at the edges of the numerical grid. Ideally, one 
would like to find a boundary condition that does not 
introduce numerical instabilities and allows gravitational 
waves to leave the grid cleanly, with no artificial reflec- 
tions. This is in itself a very difficult problem, since in 
the first place, there is no local boundary condition that 



where r = ^/ x 2 + y 2 + z 2 and where the constant 
fo is taken to be one for diagonal metric compo- 
nents and zero for everything else. The radiative 
condition assumes that the boundaries are in the 
wave zone, where the speed of light is essentially 
one, and where the gravitational waves behave as 
spherical wavefronts. This boundary condition has 
been used before by other authors [7j,|l0||, and it 
has been found that in practice it is very good at 
absorbing waves. 

It is in fact easier to implement a differential form of 
the radiative boundary condition than to use fl27|) 
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directly. Consider a boundary that corresponds to 
a coordinate plane Xi=constant. Condition (|27| ) 
then implies: 

^d t f + d i f + ^(f-f )=0. (28) 

One can now use simple finite differences to im- 
plement this last condition. In our code we have 
implemented condition ( p8| ) consistently to second 
order in both time and space. 

• Robin boundary condition: This is a different type 
of "extrapolating" boundary condition, where one 
assumes that for large r a given field behaves as: 



f(x l ) =f + k/r, 



(29) 



with k constant. This condition is clearly related 
to the radiative condition described above, but it 
contains no information about the time evolution. 
Just as we did with the radiative condition, we in 
fact implement the Robin condition in differential 
form: 



Sif + - (/-/o) = 0. 
r 



(30) 



The Robin boundary condition is usually better 
suited for solving elliptic problems than for use on 
dynamical variables. 

Most of the simulations discussed below have been per- 
formed using the radiative boundary condition (53) for 
the dynamical variables, and the Robin boundary condi- 
tion (]3(]) both for constructing the initial data and for 
solving the maximal slicing condition. Whenever a dif- 
ferent boundary condition is used, we say so explicitely. 



III. APPLICATIONS 

In this section we will apply the previous system of 
conformal trace-free equations, exploring different imple- 
mentations, in a series of numerical experiments with dif- 
ferent spacetimes. The various implementations we con- 
sider are: 

• Promoting the T's to independent variables. 

• Use the momentum constraints on the evolution 
equation for the T's. 

• Enforcing trA = 0. 

• For maximal slicing, enforcing trA" = 0. 

We will study the effects of these different implemen- 
tations using strong gravitational waves spacetimes. 

All the numerical simulations presented here are car- 
ried out with the Cactus code for numerical relativity 
co-developed in our NCSA/Potsdam/Wash U collabora- 
tion and elsewhere. 



A. Vacuum Spacetimes 

We begin our discussion of the numerical simulations 
with vacuum spacetimes in this subsection, examining 
the evolution of both strong gravitational wave and black 
hole spacetimes. In particular, we use the gravitational 
wave simulations to illustrate the effects of the various 
implementations of the CT approach. 



1. Pure Gravitational Waves 

We first turn to pure gravitational wave spacetimes. 
The low amplitude linear case has been studied, with 
a full 3D code, and published previously, (a) in both 
the standard ADM formulation and the Bona-Masso hy- 
perbolic formulation by [p4f , where no fundamental dif- 
ferences were seen in performance at that time, and 
(b) by Shibata and Nakamura Q| and Baumgarte and 
Shapiro jl0| in the CT approach as described above. The 
Baumgarte and Shapiro Jl0| work particularly showed the 
strength of the CT formulation in the linearized case. 
Here we extend the study of these systems to include 
highly dynamic, strong field regimes. The study here is 
limited to tests that show the strengths and weaknesses 
of the different formulations. A study of the physics of 
collapsing waves in full 3D numerical relativity is pre- 
sented elsewhere 

We consider here a three-metric of the form originally 
considered by Brill ^(J : 

ds 2 = * 4 [e 2 " (dp 2 + dz 2 ) + p 2 d<p 2 } = ^ds\ (31) 

where q is a free function subject to certain regularity 
and fall-off conditions. Different forms of the function 
q have been considered by different authors ]27|-|30t|, but 
most work so far has concentrated only in constructing 
and analyzing the initial data. 

As in Ref . [^| , we use a generalized form for the func- 
tion q, giving it a full 3D dependence, following [pl|-p4[: 



Q 



2 _ r 2 (1 + cp 2 cos 2 (m8)) 
ape 



(32) 



where a and c are constants, r 2 = p 2 + z 2 and m is an 
integer. In this paper we focus on the axisymmetric case, 
c = 0, for simplicity, although using a non-zero value of c 
does not affect the results we discuss below. All the runs 
discussed here where performed using an iterative Crank- 
Nicholson (ICN) scheme with 3 iterations (see appendix), 
and radiative boundary conditions. 

The first case presented is an initial configuration with 
amplitude a=4, corresponding to a strong wave, but not 
quite strong enough to collapse to a black hole. In the 
evolution of this data set the wave implodes through the 
origin, oscillates a few times, and finally disperses back to 
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FIG. 1. a) Evolution of the minimum value of the lapse for 
an axisymmetric Brill wave data set with a=4, using the stan- 
dard ADM formulation with maximal slicing. The simulation 
crashes at t=8 with a catastrophic collapse of the lapse, b) 
Evolution of the maximum value of the trace of the extrinsic 
curvature K. 

infinity leaving flat space behind, but in a non-trivial spa- 
tial coordinate system p5fl . The evolution of this space- 
time is highly non-linear, and the final configuration has 
metric components with a large dynamical range. 

In Fig. [l]a we show the evolution of the minimum value 
of the lapse over the grid for a simulation done with the 
standard ADM formulation, using maximal slicing, no 
shift and a radiative boundary condition. For this par- 
ticular simulation we used a resolution of Aa;=0.08 and 
67 3 grid points. Also, we used the fact that our data 
is symmetric across coordinate planes to evolve only one 
octant. The simulation crashes at t ~ 8 when the lapse 
collapses catastrophically in response to a blow up of the 
extrinsic curvature. Fig. |l|b shows the evolution of the 
maximum value of the trace of the extrinsic curvature K. 
Notice that even though we are using maximal slicing, K 
does not remain zero, and blows up towards the end of 
the simulation. The fact that K does not remain zero 
is not surprising, since the maximal slicing condition is 
solved numerically, and thus a residual time derivative of 
K is to be expected. The catastrophic blow-up, however, 
is a different matter and points towards the existence of 
an unstable solution of our system of equations. 

Fig. ^ shows the same simulation, but now using the 
so-called "K-driving" technique Q . The idea here is to 
add counter terms to the elliptic equation for the lapse 
to drive the numerically produced non-zero K (the trace 
of the extrinsic curvature) back towards zero. With K- 
driving, K remains much smaller until close to the point 
of crashing at t=9, with a catastrophic blow-up of the 
lapse at the end. This shows that a better control of 
the time slicing is not enough to cure the instability in 
the evolution: There exist unstable modes that are not 
controlled by keeping the value of K small. (For an anal- 
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FIG. 2. a) Evolution of the minimum value of the lapse for 
an axisymmetric Brill wave data set with a— 4, using the stan- 
dard ADM formulation with maximal slicing and a K-driver. 
The simulation goes somewhat further, and now crashes at 
t=9 with a catastrophic blow-up of the lapse, b) Evolution of 
the maximum value of the trace of the extrinsic curvature K. 
The trace now remains much smaller during the simulation. 

ysis of possible unstable modes of the ADM equations, 
see @.) 

Next, we show the evolution of the same system using 
again maximal slicing, and different implementations of 
the CT formulation. In Fig. || we show again the central 
value of the lapse for the same initial data. The different 
runs correspond to the following cases: 

use of use momentum force remove 
Y % constraints K=0 trA 



Res 


no 




no 


no 


Gam 


yes 


no 


no 


no 


Mom 


yes 


yes 


no 


no 


AFK 


yes 


yes 


yes 


no 


AFA 


yes 


yes 


no 


yes 


AF2 


yes 


yes 


yes 


yes 



The first run uses the implementation denoted "Res" 
(for rescale) . It differs from the standard ADM equations 
only in the conformal rescaling and the fact that <p and 
K (which enter into the evolution equation for Aij) are 
now evolved separately. The second run, with the im- 
plementation denoted "Gam" (for gamma), introduces 
the r 4 , but does not use the momentum constraints to 
rewrite their evolution equations. The third run uses the 
implementation "Mom" (for momentum constraints) and 
represents a straightforward coding of the the full set of 
CT equations [ffljlCfl, where the momentum constraints 
are used to modify the evolution equations for the T l , 
but without adding anything else. In the fourth run, 
which uses the implementation "AFK" (for "actively en- 
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FIG. 3. Evolution of the minimum value of the lapse for 
an axisymmetric Brill wave data set with a=4, using the 6 
different variations of the CT system described in the text. 

forcing K" ) , we have forced K to remain zero by simply 
not evolving it, and we have also kept </> time independent 
(see Eq. <fivj)). In the fifth run we use the implementation 
"AFA", where we allow K to evolve freely, but actively 
force A (the trace of Aij) to remain zero by subtracting 
it from after each time step: 



Ai 



Aij - tvA. 



(33) 



And finally, in the sixth run we use the implementation 
"AF2" that combines implementations AFK and AFA 
above by actively enforcing both K=0 and A=0. Notice 
that both K and A should be zero in principle in an exact 
evolution using the CT equations with maximal slicing, 
but they do not remain so in actual numerical evolutions 
unless actively enforced. 

As can be seen from the figure, runs Res, Gam, Mom, 
AFK and AFA eventually crash, but run AF2 with dou- 
ble active enforcement does not, at least for the time 



scale under study. The lapse returns to unity, and the 
final static spacctime can be followed for a long time with 
no sign of an instability (we have in fact followed run AF2 
past i=100 and it still remains stable). From the figure 
we also see that by enforcing only K=0 or A—O sepa- 
rately, as is done in runs AFK and AFA, one still obtains 
improved stability, with the simulations crashing at late 
times after the lapse has already returned to 1. This 
shows that by enforcing only one of the two constraints, 
and keeping the other options turned on, we still get a 
rather robust system when compared to standard ADM. 
Moreover, enforcing A—O appears to be more important 
than enforcing K—0, as can be seen from the fact that 
run AFA crashes much later than run AFK. 

Finally, notice that run Gam crashes even sooner than 
run Res, which shows that it is in fact better not to use 
the T l than to use them without modifying their evo- 
lution equation. For understanding the need to use the 
momentum constraints in the CT approach, see the com- 
panion paper Q. 

We note that the results found above for the different 
implementation are generic for strong gravitational wave 
spacetimes, quite independent of the precise parameter 
choices. However, for weak gravitational waves in the lin- 
ear regime, the straightforward coding of the CT equa- 
tions (implementation "Mom") leads also to stable evo- 
lutions as do the AFK, AFA and AF2 cases. In Fig. [| we 
show again the minimum value of the lapse for the evo- 
lution of a wave with an amplitude of a = 0.01, using the 
ADM formulation and also the Mom, AFK and AFA im- 
plementations of the CT system (since the lapse remains 
very close to 1, we are in fact plotting (a — 1) x 10 5 ). 
We see that while the ADM run crashes at an early time 
(t ~ 15) with a catastrophic collapse of the lapse, all 
three implementations of the CT equations give stable 
evolutions and yield basically the same results for a weak 
wave. We have followed these three runs past t=100 with 
no instabilities developing (the AF2 implementation is in 
fact just as stable, but we don't include it in the figure). 

From these studies (and many others with different 
parameters that we have done) we can conclude that, for 
maximal slicing, the CT formulation has better stability 
properties for the evolution of strong field systems, as 
long as: 

• The T l are promoted to independent variables. 

• The momentum constraints are used to transform 
the evolution equation for the T l . Evolving the 
without modifying their evolution equation is worse 
than not using them at all. 

• The trace of the extrinsic curvature K is actively 
forced to be zero (the definition of maximal slicing) . 

• The trace of the A^ is also actively forced to be 
zero. 
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FIG. 4. Evolution of the minimum value of the lapse for an 
axisymmetric Brill wave data set with a=0.01 for the ADM 
system and three variations of the CT system (Mom, AFK, 
AFA). Notice that since the lapse remains very close to 1, we 
are in fact plotting (a — 1) x 10 5 . 

So far we have focused on the issue of long term sta- 
bility. Now we want to compare accuracy of the CT 
and ADM formulations. We concentrate on the best im- 
plementation of the CT equations, the one we labelled 
AF2. In Fig. | we show the L2-norms of the hamiltonian 
constraint for the a=0.01 and a=4 cases discussed above, 
using the ADM (solid line) and the AF2 systems (dashed 
line). In both cases we see that for the ADM system, the 
L2-norm of the hamiltonian constraint grows more or less 
linearly for some time (this is more evident in the a=0.01 
case) until just before the crash when it blows up catas- 
trophically. In contrast, in the AF2 runs the L2-norm of 
the hamiltonian constraint initially grows faster, but it 
later settles on a constant value. The fact that the ADM 
runs are more accurate than the AF2 runs at early times 
appears to be quite generic: we have found essentially 
the same behavior for all the different parameters that 
we have studied. 

We have also performed convergence tests by running 
the same initial data with different resolutions, and we 
have found that both the ADM and AF2 evolutions are 
second order accurate. As an example of this, Fig. ^ 
shows the L2-norms of the hamiltonian constraint for 
both the ADM and the AF2 systems for two different 
resolutions: The dashed lines show the L2 norm for a 
resolution of <ix=0.16 (35 3 grid points), while the solid 
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FIG. 5. L2-norms of the hamiltonian constraint for the 
a=0.01 and a=4 cases, using the ADM (solid line) and the 
AF2 systems (dashed line). 
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FIG. 6. Convergence of the L2-norms of the hamiltonian 
constraint for the a=4 case for both the ADM and the AF2 
systems. The dashed lines show the L2 norm for a resolu- 
tion for of 0.16, while the solid lines show the L2 norm for a 
resolution for of 0.08 multiplied by a factor of four. 

lines show the L2 norm for a resolution of dx=0.08 (67 3 
grid points) multiplied by a factor of four. For second 
order convergence the solid and dashed lines should fall 
on top of each other. From the figure we see that this is 
indeed true for most of the run in both cases. For the 
ADM run, second order convergence starts to fail shortly 
before the crash. On the other hand, for the AF2 run 
we obtain slightly degraded convergence (but still better 
than first order) for times between <=5 and £=15 when 
the spacetime is very dynamic, indicating that we haven't 
quite reached the second order convergence regime for the 
resolutions considered here. 

Though in this section we have concentrated in the 
case of maximal slicing, we should mention that we have 
also performed many simulations using the generalized 
"1+log" slicings. The results are in fact very similar to 
those reported here, except for the fact that implcmcn- 
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tations AFK and AF2 can no longer be used (since K in 
non-zero for these slicing conditions). We find that for 
these algebraic slicings, implementation AFA is by far 
the best performer. 

In the following subsections, we show that the above 
results on the stability and accuracy of the ADM and CT 
systems are basically the same for systems ranging from 
black holes to spacetimes coupled to dynamical source 
fields. 

2. Black Holes 

Black holes have been the target of an intense research 
effort in recent years in numerical relativity, and have 
proved particularly difficult to handle in 3D evolutions. 
In the "standard" numerical evolution of black holes us- 
ing the ADM equations together with singularity avoid- 
ing slicings, 3D simulations generally develop instabilities 
and crash before t — 50M, where M is the mass of the 
system ]36|-|3l| . This falls far short of the time required 
to model the complete inspiral of two black holes, or even 
the head-on collision. Still, singularity avoiding slicings 
combined with the ADM equations make it possible to 
evolve through a brief part of the merger phase of two 
black holes with momenta and spins, and from this point 
of view give the most generally applicable method avail- 
able. Future cures for grid stretching are expected to 
be based on black hole excision |3j],|l| or characteristic 
slicings f4l|| . 

In the following we carry out a preliminary study of 
the CT formulation in black hole evolutions with grid 
stretching. It is inevitable that the sharp peaks that de- 
velop in the metric function due to grid stretching will 
cause the code to crash at some point in the evolution. 
We consider the evolution of the Misner data as a con- 
crete example. The 3D numerical evolution of the Mis- 
ner data in the standard ADM setting with singularity 
avoiding slicing has previously been studied using the so- 
called "G" code |37|,[3l) and its derivatives Q , developed 
by the NCSA/WashU group. Comparable results for a 
single black hole can be found in Jl8| . 

In Fig. and Fig. || we compare the results of evo- 
lutions of Misner data with the separation parameter 
fi = 2.2, corresponding to two initially well separated 
black holes, on a grid of size 130 3 with grid spacing 
0.08. The only difference in the simulations is the sys- 
tem of equations used to carry out the evolution (ADM 
vs. AF2); all computational parameters, such as param- 
eters in the ICN finite differencing scheme, grid parame- 
ters, radiative boundary conditions, and maximal slicing 
condition are the same. In Fig. fj], first panel, we show 
the radial-radial metric component along a line on the 
equatorial plane at various times for the ADM case. We 
can clearly see the familiar ever-growing peak caused by 
the grid stretching associated with singularity avoiding 
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FIG. 7. Evolution of the radial-radial metric component 
along a line on the equatorial, plane at various times for Mis- 
ner data (fi = 2.2). Plots are every 3.5M in time. The ADM 
system crashes after t = 14M, while the AF2 system remains 
stable. 
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FIG. 8. Evolution of the lapse function along a line on the 
equatorial plane at various times for Misner data (n = 2.2). 
Plots are every 3.5M in time. The ADM system crashes after 
t = 14M, while the AF2 system remains stable. 

slicings. In the first panel of Fig. || we show the lapse 
function along a line on the equatorial plane at various 
times for the ADM case, and here an instability becomes 
apparent at around t = 14M which is not yet reflected 
in the metric. This short wave length instability grows 
rapidly and causes the code to crash at t — 1AM . In the 
second panel we show the AF2 case. No metric insta- 
bility is seen until towards the end of the simulation at 
t = 2AM, although the peak appears to be deformed. At 
this time the radial metric function peak has grown to 
about two times higher than that attained in the ADM 
case. The lapse for the AF2 case in Fig. || does not show 
an instability. 

However, note that a smooth and stable evolution of 
the lapse does not mean that the computed data is still 
useful. To emphasize this point, Fig. ^| shows the same 
run as above with AF2 on a smaller grid with only 66 3 
points, but with the same grid spacing as before (so the 
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FIG. 9. Evolution of the lapse and the metric at various 
times for Misner data (/i = 2.2). Plots are every 5M in time. 
With the AF2 system the evolution remains stable even after 
the metric peak is severely deformed. 
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FIG. 10. Evolution of the L2-norm of the Hamiltonian con- 
straint for Misner data (fj, — 2.2). The ADM system crashes 
at around t — 1AM, while the AF2 remains stable. How- 
ever, the accuracy of the AF2 run degrades significantly after 
around t = 20M. 

boundaries are much closer in). While ADM crashes 
when the gradients in the metric become too severe, the 
AF2 run is able to continue with a smooth lapse even 
after the metric becomes deformed (cmp. ]l8fi where the 
evolution of the metric is not discussed). The lapse even- 
tually collapses in the whole grid, freezing the evolution 
(so one could keep running "forever" , but the evolution 
becomes meaningless). 

Next, we compare the accuracy of both simulations. 
In Fig. [lC] we show the L2 norm of the Hamiltonian con- 
straint for a grid size of 130 3 . The dashed line represents 
the ADM run, and the solid line the AF2 run. We see 
that the ADM results are more accurate than the AF2 
results until just after time t — 14M, when the insta- 
bility in the ADM evolution begins to dominate and the 
code crashes (with higher resolution this crash time can 
be delayed somewhat). Starting at around t = 20M for 



AF2, there is a spurious growth in the Hamiltonian con- 
straint that corresponds to the deformation in the metric. 
For maximal slicing one expects continuous growth of a 
smooth metric peak, but with AF2 the shoulder in the 
lapse seems to overtake the outward movement of the 
metric peak, freezing its growth in an irregular manner. 

These results for black holes with grid stretching can- 
not be compared directly to the wave runs in the pre- 
vious section because in the case of the black hole runs 
we do not approach a static final state. However, the 
CT formulation still offers some advantages over ADM 
in achievable run time. We find stability far beyond were 
the runs are meaningful, and it remains to be explored 
how far one can push the CT runs while maintaining 
convergence. 

B. Matter Spacetimes 

In the previous sections we studied the stability prop- 
erties of the vacuum Einstein equations. What will hap- 
pen if these equations are coupled to dynamical matter 
sources that are themselves governed by evolution equa- 
tions coupled to the spacetime geometry? The complete 
set of equations can now have more complicated types of 
unstable modes. What would be the effects of switching 
from the ADM formulation to the CT formulation? 

To respond to this question we consider next the fol- 
lowing systems: (i) the evolution of boson stars gov- 
erned by the scalar field Klein-Gordon equation and (ii) 
the evolution of neutron stars governed by the hydrody- 
namical equations (general relativistic Euler equations). 
The numerical evolution of the Klein-Gordon equation 
is straightforward with many well-known stable schemes. 
However, the numerical evolution of the hydrodynamical 
equations is considerably more challenging, especially in 
the presence of shocks or highly relativistic flows. For 
this purpose we use a recently developed hydrodynamical 
code p2[ which employs a conservative formulation of the 
equations together with high-resolution shock-capturing 
(HRSC) schemes based on approximate Riemann solvers. 
In |l2| we demonstrated that this code is capable of han- 
dling hydrodynamical evolutions in a stable and accurate 
fashion for a range of scenarios. 

We focus here on analyzing the stability and accuracy 
of evolutions of both static boson stars and static neu- 
tron stars using the ADM formulation and the AFA im- 
plementation of the CT equations discussed above. We 
use the AFA implementation rather than AF2 because 
the simulations discussed here have all been performed 
using algebraic slicings and implementation AF2 applies 
only to maximal slicing. The main motivation for this 
has been the fact that, as we will show below, implemen- 
tation AFA with algebraic slicings already gives excellent 
results when compared with standard ADM and is far less 
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computationally expensive than runs that use maximal 
slicing. 



1. Boson Stars 

We begin with a simple kind of matter source: self- 
gravitating scalar fields. This system has served as a 
useful testbed for developing numerical techniques for 
dealing with rclativistic matter coupled to the Einstein 
equations |4^,|35|,|4|,^5| , and also has a distinguished his- 
tory in the field, having provided the first example of 
critical phenomena in relativity ]4q| . 

The dynamics of a massive scalar field are described 
by the minimally coupled Klein-Gordon (KG) equation 



U g (j> = m <f>, 



(34) 



(see, e.g. Q). The KG equation can be obtained from 
the Lagrangian 



which leads to the stress-energy tensor 

2 5C 



T, 



[lis 



-9 Sgf" ' 



(35) 



(36) 



which is used as the matter source for the Einstein equa- 
tions. 

Self-gravitating massive scalar fields have bound, star- 
like solutions called boson stars with stability properties 
very much like those of neutron stars. These objects have 
been studied numerically, extensively in ID fli^pif ] and 
also in 3D |3jJ. Apart from the fact that their evolution 
equation is much simpler than the hydrodynamical equa- 
tions, boson stars are also easier to handle numerically 
when compared to neutron stars because they have no 
sharp changes in the density distribution near the sur- 
face layer of the star. For more details on the properties 
of boson stars and their behavior under perturbations 
see Q and references cited therein. 

We perform our numerical evolutions of boson stars by 
writing the KG equation as a flux-conservative system of 
the form 



ii a = d b F b a + S b a u b 



(37) 



where u contains the scalar field and its time and space 
derivatives. The method used to integrate this equation 
is a symmetrized MacCormack with both directional and 
Strang splitting. Symmetrized here means that the or- 
der of left-hand and right-hand differencing changes every 
time step (this improves the stability of the scalar field 
evolution). The code for solving the KG equation con- 
verges to second order in time and space. See Ref. [p7|j48| 
for details of the numerical methods. 



We have carried out evolutions of equilibrium boson 
star configurations with the metric background held fixed 
artificially (not updating the metric functions), and evo- 
lutions of the metric of such configurations with the 
scalar field held fixed artificially (not updating the scalar 
field), for a range of compactness of the boson stars, us- 
ing both the ADM and AFA schemes. For all these cases, 
we have seen that the simulations are stable and second 
order convergent. The case of coupled spacetime-scalar 
field evolution is much more challenging, and we focus 
on that below. 

We begin by showing an equilibrium boson star 
with a central density near the maximum stable value 
(field strength at center cf>o = 0.26, total mass 
M = 0.6322 TO p 2 /m, with m p the Planck mass, m the 
mass of the scalar field). In Fig. [n], we show the evo- 
lution of radial metric component g rr in a fully coupled 
simulation, using a three step ICN scheme, 1+log slicing 
with N = 2, no shift, a radiative boundary condition on 
the metric, and a flat boundary condition on the scalar 
field. A 32 3 grid is used to cover only one octant. In the 
first panel we show the results of the ADM evolution. We 
see that for a short time, the spacetime remains nearly 
static (as it should). However, a short wavelength insta- 
bility becomes significant by time t=7, and quickly grows 
causing the code to crash. The time t here is expressed in 
terms of the intrinsic oscillation time scale of the scalar 
field (the exact equilibrium boson star field has the form 
ip(r)e lt ). In the second panel we show the evolution with 
exactly the same setup but using now implementation 
AFA instead of ADM. We see that the static configu- 
ration is maintained for a much longer time. Towards 
the end of the evolution, near t — 150, we see that nu- 
merical error starts to build-up near the boundary of the 
computational domain. 

In Fig. [l2] below, we compare the L2-norm of the hamil- 
tonian constraint for the ADM and AFA runs. We see 
that at early times the ADM run gives a more accurate 
result, but instabilities cause the L2-norm to blow by 
t ~ 8. For the AFA run the constraint violation is larger 
at first, but the evolution remains stable or a much longer 
time. The oscillation of the hamiltonian constraint we see 
here can be understood as a reaction of the scalar field to 
the numerical truncation error, which can be interpreted 
as a kind of perturbation. The frequency of these oscil- 
lations coincides with the ones obtained in ID studies of 
perturbed boson stars. Notice that with the ADM run 
the code crashes so early that one can not even see the 
first oscillation. 



2. Static Neutron Stars 

We turn now to the study of hydrodynamical evolu- 
tions of neutron stars. In j|2) we developed a three- 
dimensional, fully relativistic code to integrate the hy- 
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FIG. 11. Evolution ol the radial metric function g rr us- 
ing ADM (upper panel) and the AFA implementation (lower 
panel). The ADM evolution crashes at t ~ 8. 




FIG. 12. Evolution of the L2-norm of the Hamiltonian con- 
straint for a static and stable fully coupled boson star using 
the ADM and AFA systems. The evolutions where carried 
out on a 32 3 grid, with a resolution of Ax = 0.45. 



drodynamical equations coupled to the ADM equa- 
tions. Convergence studies using polytropic neutron stars 
showed that the code is second order accurate in both 
space and time. For the integration of the hydrody- 
namical equations we used HRSC schemes of the total- 
variation-diminishing (TVD) class, with a piecewise- 
linear reconstruction of a sufficient set of hydrodynamical 
variables (rest-mass density, three-velocity and internal 
energy density). For more details on the schemes avail- 
able in the code, see In the studies reported in this 
paper we use the ICN scheme for the integration of the 
spacetime equations (either ADM or AFA) and Roe's ap- 
proximate Riemann solver for the hydrodynamical equa- 
tions. We use "1+log" slicing with N = 2. 

As in the boson star studies we have first considered 
evolutions which test separately the individual compo- 
nents of the code. In these, we either solve the hydrody- 
namical equations in a prescribed (static) spacetime or 
the gravitational field equations for a prescribed matter 
source. In particular, we have evolved static neutron star 
configurations with a zero-temperature polytropic equa- 
tion of state, of the form P = Kp T (where P is pressure 
and p is rest- mass density). This included stars with a 
large polytropic index T (very stiff) having density pro- 
files with a discontinuous first derivative at the surface. 
In the case of prescribed matter sources, we have con- 
firmed that the comparison of the AFA and AF2 systems 
to the ADM system, in terms of stability and accuracy, 
remains the same as in the vacuum cases studied above. 
Static neutron stars with polytropic index T = 2 have 
also been studied in |ll| using the CT equations with 
prescribed hydrodynamical sources. 

We focus next on the coupled spacetime and hy- 
drodynamical evolution of static Tolman-Oppcnheimer- 
Volkoff (TOV) |^9| neutron stars (in isotropic coordi- 
nates). Again, we compare the results obtained using 
the AFA implementation of the CT equations to those of 
the ADM equations. In principle both the matter distri- 
bution inside the star and the spacetime should remain 
static. In practice they evolve due to truncation errors 
of the finite-difference scheme, with the hydrodynamics 
and the spacetime responding to one another. The static 
TOV solution provides a reference to monitor the accu- 
racy of the coupled numerical evolution. Note that in 
these evolutions, static outer boundary conditions were 
used. 

In Fig. [l^, we show the evolution of the L2-norm of the 
Hamiltonian constraint for a polytropic, N = 1, TOV 
star of gravitational mass 1.4M Q and compactness ratio 
M/R = 0.146. A 64 3 grid is used to cover the first oc- 
tant, with dx — dy = dz = 0.34km. The dashed line 
corresponds to the ADM system and the solid line to the 
AFA system. Again, as in the vacuum studies, we see 
that the ADM evolution suddenly becomes unstable at 
roughly 2.7ms, while the AFA evolution remains stable 
after more than 6ms (we followed the evolution for more 
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FIG. 13. Evolution of the L2-norm of the Hamiltonian con- 
straint for a N=1.0 polytropic neutron star model (coupled 
spacetime and hydrodynamical evolution) . The ADM system 
crashes after less than 2.7ms, while the AFA system evolves 
stably for a significantly longer time. A 64 3 -grid was used to 
cover the first octant. 



than twice that). 

In Fig. [l4| we show the evolution of the radial compo- 
nent of the metric (constructed from the evolved Carte- 
sian metric components). The first panel of Fig. |lj cor- 
responds to the evolution obtained with ADM. We see 
that the star basically maintains its initial equilibrium, 
until the high-frequency instability crashes the code. In 
the second panel, we show g rr at various times, obtained 
with the AFA implementation. All other parameters are 
the same as in the ADM evolution. The ADM run is more 
accurate, before it becomes unstable, while the AFA run 
is stable but less accurate (there is a secular drift away 
from the initial configuration). 

The truncation errors of the coupled evolution code ini- 
tiate a pulsation of the star in, mainly, its radial modes 
of pulsation. These pulsations are damped in time due to 

The 



the viscosity of the numerical scheme (see [50 51 
TVD schemes we are using describe well the physical 
pulsations of the fluid, except in a small region around 
the center of the star, where short wavelength noise ap- 
pears in the radial velocity. Our trials with other HRSC 
schemes show that this behavior seems to be generic for 
higher order HRSC schemes In all such schemes, the 
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FIG. 14. Comparison of the evolution of the radial metric 
component for a N = 1.0 polytrope with ADM (left panel) 
and AFA. The evolution with the latter system proceeds well 
beyond the time at which the ADM system becomes unstable. 

radial momentum near the center has a small residual 
value of constant sign. This momentum appears in the 
r.h.s. of the evolution equation for f 1 (Eq. (po|)). This, 
in turn, leads to an error in the spacetime evolution. It 
is noteworthy that this does not cause an instability in 
the coupled evolution, except at very late times, when 
the violation of the Hamiltonian constraint has already 
become extremely large. 

We note that as the TVD schemes are only first- 
order accurate at local extrema, such as the maximum 
of the density at the center of the star, so the increase 
in the Hamiltonian constraint at the center converges 
to roughly first order with increasing resolution. Away 
from the center, the scheme is second order convergent. 
The convergence of the L2-norm of the Hamiltonian con- 
straint with the AFA system, for different grid-sizes (and 
for the same initial configuration as above), is shown in 
Fig. P. 



IV. DISCUSSION AND CONCLUSIONS 

In this paper we have studied the stability of three- 
dimensional numerical evolutions of the Einstein equa- 
tions in a formulation that separates out the conformal 
and traceless parts of the system. In our study we have 



*We have extensively experimented with other hydrodynam- 
ical evolution schemes. If one uses a first-order (Godunov) 
scheme, using piecewise constant reconstructed data for the 
Riemann problem, instead of piecewise linear, the radial ve- 
locity oscillates around zero near the center of the star, with- 
out any short wave length noise. But, a low-order scheme is 
not capable of accurately describing the evolution of the stel- 
lar surface where the density distribution is changing rapidly 



(unless prohibitively large grids are used) and large errors 
from the surface layers soon propagate inside the star. We 
have also experimented with a mixed system: first-order near 
the center and second-order near the surface. In this case the 
error grows at the interface of the two regime, yielding a even 
less accurate evolution overall. 



13 



1.0 



0.8 



„ E 0.6 - 



E 0.4 



0.2 



0.0 



16 

32 

64 



0.0 



1.0 



2.0 



3.0 

t (ms) 



4.0 



5.0 



6.0 



FIG. 15. Convergence of the L2-norm of the Hamiltonian 
constraint, at three different resolutions, for a N=1.0 poly- 
tropic neutron star model. The AFA system is used. 



considered different spacetimes including gravitational 
waves, black holes, boson stars and neutron stars. 

We investigated several implementations of the 
conformal-traceless (CT) evolution equations. We iden- 
tified two of them which give the best long term stabil- 
ity behavior: the AF2 implementation for maximal slic- 
ing, and the AFA for algebraic slicings. The AFA imple- 
mentation actively enforces the trace of the conformally 
rescaled extrinsic curvature (A) to zero at each step of the 
time evolution, while the AF2 implementation enforces as 
well the fact that the trace of the extrinsic curvature (K) 
should vanish in maximal slicing. On the analytic level, 
the CT evolution equations imply that A = throughout 
the evolution, but this is inevitably violated in numeri- 
cally evolution due to truncation error, unless actively 
enforced. Similarly, for maximal slicing, K will not re- 
main zero numerically unless actively forced to do so. 
We find that these two implementations of the CT equa- 
tions lead to a more stable evolution compared to what 
one can obtain using the standard ADM evolution equa- 
tions, under the same resolution, boundary condition and 
grid parameter choices, for all systems investigated. In 
comparison, a straightforward implementation of the CT 
equations ("Mom") is capable of giving a stable evolu- 
tion for weak but not strong field systems. We should 
also mention that we have recently become aware of the 
work of Lehner, Huq and Garrison [p~5[ where a compar- 
ison of the ADM and CT formulations has been carried 
out and where it is also found that freezing the evolution 
of K (what these authors call "locked evolution") im- 
proves considerably the stability of simulations that use 
the CT formulation. 

Beyond stability, we have also compared the accuracy 
of the evolutions obtained by the ADM equations and 
CT equations. For all spacetimes considered we have 
found that the ADM system is consistently more accu- 



rate than the CT system in short term evolutions, before 
the instabilities set in. Although at present we can offer 
no explanation of this difference in accuracy between the 
different formulations, we believe that it is not a conse- 
quence of our numerical implementation, but is rather a 
property of the system of differential equations. It there- 
fore points in the direction for a possible improvement of 
the CT approach. We note that formulations combining 
the CT approach and the hyperbolic approach have been 
proposed |5^,|5J|. A similar investigation of the stabil- 
ity and accuracy properties of such formulations will be 
presented elsewhere. 

In this paper we have focused on the implementations 
and the numerical properties of their evolutions. Some 
understanding of the different stability of properties on 
the analytic level is discussed in a companion paper B. 
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APPENDIX A: STABILITY ANALYSIS OF THE 
ITERATIVE CRANK-NICHOLSON SCHEME 

The numerical scheme used for the simulations de- 
scribed in this paper is the so-called iterative Crank- 
Nicholson (ICN) scheme, which is an iterative, explicit 
version of the standard implicit Crank-Nicholson (CN) 
scheme p4] , p5| . The idea behind this method is to solve 
the implicit equations by an iterative procedure, where 
each iteration is an explicit operation depending only 
on previously computed data. Normally, this process 
is stopped after a certain number of iterations, or until 
some tolerance is achieved. For a linear equation (and in 
particular in one dimension), the iterative procedure can 
easily be much more computationally expensive than the 
matrix inversion required to solve the original implicit 
scheme. For a non-linear system, however, solving the 
implicit scheme directly can prove to be extremely diffi- 
cult. 



14 



In this appendix we study the stability properties of 
the ICN scheme in the particular case of the simple wave 
equation, and derive two very important results: 

• In order to obtain a stable scheme one must do 
at least three iterations, and not just the two one 
would normally expect (two iterations are enough 
to achieve second order accuracy, but they are un- 
stable!). 

• The iterative scheme itself is only convergent if the 
standard Courant-Friedrichs-Lewy (CFL) stability 
condition is satisfied, otherwise the iterations di- 
verge. 

These two results taken together imply that there is 
no reason (at least from the point of view of stability) to 
ever do more that three ICN iterations. Three iterations 
are already second order accurate, and provide us with a 
(conditionally) stable scheme. Increasing the number of 
iterations will not improve the stability properties of the 
scheme any further. In particular, we will never achieve 
the unconditional stability properties of the full implicit 
CN scheme, since if we violate the CFL condition the 
iterations will diverge. [j] 

For our stability analysis we will consider the simple 
wave equation in N-dimensions. Numerical experiments 
have shown that the full Einstein equations have essen- 
tially the same stability properties. 

Consider then the N-dimcnsional wave equation writ- 
ten in "3+1 like" form: 



8t<f> = A, d t A = J2 9 i 



(Al) 



i=i 



For the finite difference approximation to these equa- 
tions we employ the usual notation 



:=f(xi =miAx,t = nAt), 



(A2) 



with n and m = (mi, m^r) integers. The implicit CN 
scheme is then given by 



i.n+1 _ in 



At 



A 



n+l 



A 7 ' 



At 



n+1 + AH 



N 



2(Ax) 2 ^ 



(A3) 



(A4) 



tAs we were finishing this manuscript we became aware of a 
paper by S. Teukolsky were he does essentially the same anal- 
ysis and obtains the same results |Q. His analysis and ours 
complement each other, since he considers any finite number 
of iterations, while we consider only 1, 2 and 3 iterations. 
On the other hand, here we also consider the question of the 
convergence properties of an infinite number of iterations. 



where the finite difference operators 5f are defined as 

$i /mi := fnii + 1 ~ 2fmi + /m ; -l- (A5) 



The implicit CN scheme is well known to be uncondi- 
tionally stable for the wave equation (i.e. stable for any 
value of At) . 

The ICN scheme is defined in the following way 



... , 2, + Ai Al 

N 

A^=A^ + AtJ2^, 



N 



i=l 



^ (aw 

2 V m 



AH 



(A6) 
(A7) 

(A8) 



A(i) — A n i 

rn in > 



At 



2(Ax) 2 



N 

2 E 



5f + d) , (A9) 



and finally, 



6 n+1 = 

rm I'm ) 



(A10) 
(All) 



From these expressions it is clear that if the iterations 
converge, we will recover the implicit CN scheme. 

For the stability analysis of the ICN scheme we use the 
standard von Neumann ansatz |5lj57|l 



rn t \n i(k.-m)Ax 

= ^ 2 A n e i(k ' m)Ax , 



(A12) 
(A13) 



with k the "wave vector" . Notice that the highest wave 
number that can be represented on the finite difference 
grid corresponds to kiAx = tt. The stability condition 
for our numerical scheme will then be 



|A|<1- 



(A14) 



Let us consider first the "1-step" ICN scheme, that 
is, the so-called forward-time centered-space (FTCS) 
scheme. This scheme is well known to be only first order 
accurate, and unconditionally unstable. The fact that is 
only first order accurate can be easily seen from a sim- 
ple Taylor expansion in time. For the st ability analysis 
we substitute the von Neumann ansatz ( A13 ) into the 
ICN scheme defined above with « max = 1. Doing this we 
obtain 



A 2 - 2A + 1 + 2p 2 u 2 = 0, 
where p := At/ Ax is the Courant parameter and 



JV 

i=i 

u 2 := 1 — cos{kiAx) 



(A15) 

(A16) 
(A17) 



Solving for A we find 
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X = l±iV2 



pu, 



which implies 



(A18) A = 1 - p 2 u 2 ± i V2 pu 1 1 - p V/2 1 

which now implies 



(A28) 



|A| = 1 



2p 2 u 2 > 1. 



(A19) 



Comparing with ( A14 ) we conclude that the 1-step 
scheme is unstable for any value of At. 

Let us now consider the 2-step scheme. If we take the 
ICN scheme above with i max = 2, and do the appropriate 
substitutions we find 



in+l 



At A" 



P_ 

2 



N 



i=l 



A 



n+1 



2 N 
i=l 



AtA" 



(A20) 
(A21) 



As before, a simple Taylor expansion shows that this ap- 
proximation is now second order both in time and space. 
Using again the ansatz ( |A13 ) we find now that 



A 2 + 2A (p 



2„,2 



4 4 
p U 



Solving again for A we obtain 



A = 1 



2 2 
p U 



±iV2 



pu, 



which implies 



|A| 



i + pV > i. 



(A22) 



(A23) 



(A24) 



Comparing again with ( A14 ) we conclude that the 2- 
step ICN scheme is also unstable for any value of At. 
This result is surprising, since a priori one might expect 
that the 2-step scheme should behave like a predictor- 
corrector scheme, and should therefore be stable. 

Finally, let us consider the 3-step scheme. By taking 
the ICN scheme above with i max = 3, and doing the 
appropriate substitutions we now find 



N 



frn+1 = An 
m rm 



An+1 _ An 

Sim In 



+ AtA n m + "-Y,5U^ n m + AtAl l ), (A25) 



£1 
2 



+ 



AAx 



f N \ 2 

Em 



(A26) 



A Taylor expansion now shows that this 3-step scheme is 
still only second o rder accurate in both time and space. 
Using the ansatz ( A13 ) on this scheme we now find 



-, 4 4 
1 — p U 



A 2 + 2A (p 2 u 2 - 1 
And solving for A we obtain 



1 



-A 6 = o. 



(A27) 



|A| = l-pV + ^A 6 . 



(A29) 



Comparing now with (A14) we obtain the following 
stability condition 



p 2 u 2 < 2. 



(A30) 



And finally, from the fact that the maximum value of u 2 
is 2y/~N we find 



(A31) 



p < l/VN 



Notice that this is just the standard CFL condition in N 
dimensions. We then conclude that in order to obtain a 
(conditionally) stable scheme we need to do at least three 
iterations. 

Next, we address the question of the stability of the 
iterations themselves, that is, if we iterate an infinite 
number of times do we converge to something (that is, 
to the implicit CN scheme)? For this we consider two 
consecutive iteration steps (i — l,i), and subtract them 
to get 







_ 4 — !) = 



2 \ m in 

N 

M T (V*- 13 - 

2(Aa;) 2 ^ V^ m 



(i-2) 



(A32) 
(A33) 



Let us now define F±Q :— 4>m — 13 and 
i^m := — ^m" 1 '- The above equations become 



At 



N 



2{Ax) 



We now use the von Neumann ansatz again 



£ \ % „i(k-m) Ax 
J 2 A e ' , 



(A34) 
(A35) 

(A36) 
(A37) 



Substituting this ansatz back into the equations above 
we find 



A 



1 



2 2 
p U 



0, 



from which we obtain 



A = ±i—=. 
\/2 



(A38) 



(A39) 



In this case, the condition for the iterations to converge 
implies that the norm of the successive differences should 



1G 



go to zero, which in turn implies |A| < 1. Using again the 
fact that the maximum value of u 2 is 2v / 7V we see that 
the convergence condition reduces to 



p< 1/y/N. 



(A40) 



This is again the standard CFL stability condition. So 
we have just shown that if this condition is violated, the 
iterations will fail to converge. This means that there is 
no reason to try to iterate to convergence in the hope of 
improving stability. If At was too big in the first place 
the iterations will never converge. 
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